
rm(list=ls())
setwd("C:/Users/torwig/Dropbox/Samarbeid/T&Sdelemappe/Faglig/GenderConflict/Data")
library(foreign)
library(devtools)
library(margins)
library(dummies)
library(MASS)

#####################################################################################
################              FIGURE 2 - MARGINAL EFFECT            #################
#####################################################################################
mydata <- read.dta("forsimulation.dta")




#year dummies

#baseline model (region-dummies)
summary(mod1 <- glm(conflict_onset ~ v2x_gender + v2x_polyarchy
                    + v2x_polyarchy_sq + log_gdpPC + log_pop +  peaceyears + peaceyears2 + peaceyears3 + as.factor(year)-1
                  
                    
                   ,data=mydata,
                    family=binomial(link="logit")))



beta <- coef(mod1)
beta <- beta[1:8]
covvar <- vcov(mod1)[1:8,1:8] #the relevant covariance matrix (without the year fixed effects)


precolseq <- seq(from=0, to=1, by=.10)
m <- with(mydata, cbind(precolseq,
                      mean(v2x_polyarchy, na.rm=T),
                      mean(v2x_polyarchy_sq, na.rm=T),
                      mean(log_gdpPC,na.rm=T),
                      mean(log_pop, na.rm=T),
                      mean(peaceyears, na.rm=T),
                      mean(peaceyears2, na.rm=T),
                      mean(peaceyears3, na.rm=T)
                   
                      
                      
                      ))


#simulating 
beta.sim <- mvrnorm(10000, beta,covvar)
xb <- m %*% t(beta.sim)
prob <- 1/(1+exp(-xb))
predprobs <- t(apply(prob,1,quantile, probs = c(.10,.5,.90)))


pdf("simulated1.pdf", width=11,height=10)
par(mfrow=c(1,1))
plot(x=c(0,1),y=c(0,0.6),type="n",xlab="Female Empowerment index",ylab="Probability of civil conflict onset",main="Simulated probability of civil conflict onset", axes=FALSE)
axis(side=1)
axis(side=2)
j <- length(predprobs)

lengths <- seq(from=0, to=1, by=.10)

polygon(c(lengths
          , rev(lengths))
        , c(predprobs[,1]
            , rev(predprobs[,3])),
        col="orange", border=NA)
lines(precolseq,predprobs[,2],lty=1, col="black")
dev.off()
